Oregon’s marine dissolved oxygen standard states:

  1. For ocean waters, no measurable reduction in dissolved oxygen concentration may be
    allowed.

Currently, no marine DO assessment methodology exists. DEQ is in the process of planning and assembling a technical work group to help develop assessment methodology for marine OAH issues.

There are two data sets that the assessment team analyzed. 1. Continuous DO data collected at marine mooring stations (OOI data) and 2. Discrete DO samples collected on cruises along the Newport line. The following in a summary of the analysis performed on the data.

OOI data

Download continuous DO data from AWQMS:

ooi_DO_data <- AWQMS_Data_Cont(org = "OOI_(NOSTORETID)",
                          char = "Dissolved oxygen (DO)")

This data represents one station at CE01ISSM for date range 2018-04-03 - 2020-12-31.

The first step the assessment team took was to visualize the data. Chan et. al. 2019 defines some conditions for marine DO.

Status Value units
Hypoxic 1.99864 mg/l
Severe Hypoxia 0.71380 mg/l
Suboxic 0.14276 mg/l
Anoxic 0.00000 mg/l
ooi_DO_data_graph <- ooi_DO_data %>%
  mutate(result_conv = Result_Numeric* (1.4/63.9)*1.42905,
         result_unit = 'ml/l',
         date = ymd(Result_Date),
         month = month(date),  
         yrmon = as.yearmon(date), 
         year = year(date))



q <- 
  ggplot(data = ooi_DO_data_graph) +
    geom_point(aes(x = date, y = result_conv, color = Depth)) + 
    labs(title = "Dissolved Oxygen",
         subtitle = paste0(ooi_DO_data_graph$MLocID[1], ": ", ooi_DO_data_graph$StationDes[1]),
         caption = paste0("Lat/Long: ",
                          ooi_DO_data_graph$Lat_DD[1], ", ",ooi_DO_data_graph$Long_DD[1] ),
         x = element_blank(),
         y = "mg/L")+
    theme_bw()+
  geom_line(data = water_status, 
            aes(x = Date, y = Value, linetype = Status)) +
  scale_color_brewer(palette="Paired")

q

There appears to be some seasonality to the DO pattern, with regular drops into Severe Hypoxia conditions and occasional drops into Suboxic conditions.

Next, we wanted to see if there was any trend that could be seen in the data. Since there is some obvious seasonality to the data set, the assessment team ran a Seasonal Mann Kendall test on the monthly average DO concentration, using months as the seasonal component. To ensure data completeness we only included months with less than 3 missing days in the data set.

month_average <- ooi_DO_data %>%
  mutate(result_conv = Result_Numeric* (1.4/63.9)*1.42905, #convert to mg/L
         result_unit = 'mg/l',
         date = ymd(Result_Date),
         month = month(date),  
         yrmon = as.yearmon(date), 
         year = year(date),
         loc_depth = paste0(MLocID, "-", Depth,Depth_Unit )) %>%
  group_by(MLocID,loc_depth, year, month, Depth) %>%
  summarise(mon_avg = mean(result_conv),
            num_days = n_distinct(date),
            count = n()) %>%
  ungroup() %>%
  complete(loc_depth, year, month) 


months_exclude <- month_average %>%
  mutate(mo_days = days_in_month(month),
         diff = mo_days - num_days,
         exclude = ifelse(diff > 3 | is.na(mon_avg), "yes", "no")) %>%
  select(MLocID, loc_depth, year, month, exclude)


month_average_assess <- month_average %>%
  left_join(months_exclude, by = c("loc_depth", "year", "month", "MLocID"))%>%
  filter(exclude == "no") %>%
  select(MLocID, loc_depth,  Depth, year, month, mon_avg) %>%
  mutate(date = paste0(year, "-", month, "-01")) %>%
  mutate(date = ymd(date)) %>%
  mutate(variable = as.factor("DO")) %>%
  select(date, MLocID,loc_depth,  Depth, variable, mon_avg) %>%
  rename(value = mon_avg,
         site = MLocID) 



#Loop through different depths 

#create empty list to accept test results
kendall_list <- list()


for(i in 1:length(unique(month_average_assess$loc_depth))){
  
  #select unique site ID and arrange in chronological order
  #fill in NAs for missing months
  
  mo_avg_seamannkenn <- month_average_assess %>%
    filter(loc_depth == unique(month_average_assess$loc_depth)[i]) %>%
    mutate(Depth = as.numeric(Depth))
  
  #create wqdata format table
  wqdataframe <- wqData(mo_avg_seamannkenn, c(1, 2, 4), 5:6, site.order = TRUE, type = "long",
                        time.format = "%y-%m-%d")
  
  #create time series
  timeseries <- tsMake(wqdataframe, focus= wqdataframe$site[1])
  
  res <- seaKen(timeseries)
  
  #save results as a dataframe and add MLocID and n
  df <- data.frame(as.list(unclass(res))) %>%
    mutate(loc_depth = unique(month_average_assess$loc_depth)[i],
           n = length(timeseries))
  
  
  
  
  #bind results list
  kendall_list[[i]] <- df
  
  
  
} #end of for loop

#bind results of seasonal mann kendall test to single dataframe
kendall_results <- bind_rows(kendall_list) 


#trend is detected when the tau > critical value and p. value (sl) is <0.10
#reorder and rename rows for better exporting
kendall_results <- kendall_results %>%
  mutate(significance = ifelse(p.value < 0.10 & sen.slope < 0, "Signifigant (-)", 
                               ifelse(p.value < 0.10 & sen.slope > 0, "Signifigant (+)",
                                      "No Trend"))) %>%
  select(loc_depth, significance, p.value, sen.slope)


kendall_results %>%
  kbl() %>%
  kable_classic(full_width = F, html_font = "Cambria")
loc_depth significance p.value sen.slope
CE01ISSM-25.00m No Trend 0.6674365 -0.0210354
CE01ISSM-7.00m No Trend 0.5790997 0.1116678

In this data set, there is no significant trend detected. However, this is a limited data set, with only a few years covered. This data represents what was submitted to us. In future assessments, if we have a longer data set, it is possible that a trend may emerge. To visualize the data:

station <- ooi_DO_data %>%
  select(OrganizationID, MLocID, StationDes,Lat_DD, Long_DD, AU_ID ) %>%
  distinct()


#point graph
mon_avg_graph <- month_average %>%
  mutate(yearmon = as.Date(paste0(year,"-",month,"-1"))) %>%
  mutate(moname = month.name[month]) %>%
  filter(!is.na(mon_avg)) %>%
  left_join(station, by = "MLocID")

#join trend analysis to data
mon_avg_graph_trend <- mon_avg_graph %>%
  left_join(kendall_results, by = "loc_depth") %>%
  left_join(months_exclude, by = c("loc_depth", "year", "month", "MLocID")) %>%
  filter(exclude == "no") %>%
  select(-exclude)



p <- ggplot(data = mon_avg_graph_trend)+
  geom_point(aes(x = yearmon, y = mon_avg, color = Depth), size = 2, position = position_dodge(0.15)) +
  geom_line(aes(x = yearmon, y = mon_avg, color = Depth)) +
  labs(title = "Average Monthly Dissolved Oxygen",
       subtitle = paste0(mon_avg_graph$MLocID[1], ": ",mon_avg_graph$StationDes[1] ),
       caption = paste0("Lat/Long: ",
                        mon_avg_graph$Lat_DD[1], ", ",mon_avg_graph$Long_DD[1] ),
       x = element_blank(),
       y = "mg/L")+
  theme_bw()

#graph no trend values
if(mon_avg_graph_trend$significance[1] == "No Trend"){
  p = p + annotate("text", label = "No Trend", 
                   x =  as.Date('2020-06-01'),
                   y = 1.75,
                   colour = "black", size = 3.5) 
  
} #end of no trend if statement

if(mon_avg_graph_trend$significance[1] != "No Trend"){
  p = p +  annotate("text", label = "Significant Trend (p-value < 0.10)", 
                    x = ymd('2020-06-01'),
                    y = 1.7,
                    colour = "black", size = 3.5) 
  
  
  
  
  #plot trend line. This is taken from the coho trends process
  slope <- kendall_results[kendall_results$MLocID == unique(sdadm_raw_trend$MLocID)[1], "sen.slope"]
  x.delta <- as.numeric((max(mon_avg_graph$year) - min(mon_avg_graph$year)))/2
  SK.min <- mean(mon_avg_graph$sdadm, na.rm = TRUE) - x.delta*slope
  SK.max <- mean(mon_avg_graph$sdadm, na.rm = TRUE) + x.delta*slope
  
  p <- p + geom_segment(aes(x = ymd('2018-01-01'), y = SK.min,
                            xend =ymd('2021-02-01'), yend = SK.max, linetype = "SeasMannKendall Trend"),
                        size = 1.05, color = "gray49") +
    scale_linetype_manual(values=c("dotted")) +
    guides(linetype=guide_legend(title = element_blank(), order = 2))
  
  
}



#Add hypoxia lines

water_status <- data.frame(
  stringsAsFactors = TRUE,
  Status = c("Hypoxic","Hypoxic",
             "Severe Hypoxia","Severe Hypoxia","Suboxic",
             "Suboxic","Anoxic","Anoxic"),
  Date = c("1/1/2018","2/1/2021","1/1/2018",
           "2/1/2021","1/1/2018","2/1/2021",
           "1/1/2018","2/1/2021"),
  Value = c(1.4, 1.4, 0.5, 0.5, 0.1, 0.1, 0, 0)
) %>%
  mutate(Date = mdy(Date),
         Value = Value*1.4276) #convert to mg/L

water_status$Status <- factor(water_status$Status, levels = c("Hypoxic", "Severe Hypoxia", "Suboxic", "Anoxic"))    

p2 <- p +
  geom_line(data = water_status, 
            aes(x = Date, y = Value, linetype = Status)) +
  scale_color_brewer(palette="Paired")

p2

Due to the regular drops into Severe Hypoxia conditions and occasional drops into Suboxic conditions, combined with the lack of a clear downward trend, the data supports this marine assessment unit remaining as category 3B- insufficient data; potential concern.

Newport Line Data

The second data set we analyzed is ship based discrete DO samples taken along the Newport line. We loaded all the Oregon Territorial Waters (3 miles out) data we could find into AWQMS.

newport_line_awqms <- AWQMS_Data(org = "NOAANEWPORTLINE_(NOSTORETID)",
                                 char = "Dissolved oxygen (DO)")

This is a much longer dataset and represents 2 monitoring locations beginning in 1998at a number of depths.

newport_line_summary <- newport_line_awqms %>%
  mutate(SampleStartDate = ymd(SampleStartDate),
         act_depth_height = as.numeric(act_depth_height)) %>%
  group_by(MLocID) %>%
  summarise(min_date = min(SampleStartDate),
            max_date = max(SampleStartDate),
            min_depth = min(act_depth_height),
            max_depth = max(act_depth_height),
            num_depths = n_distinct(act_depth_height)
            )


newport_line_summary %>%
  kbl() %>%
  kable_classic(full_width = F, html_font = "Cambria")
MLocID min_date max_date min_depth max_depth num_depths
NH01 1998-08-14 2021-03-24 0 35 96
NH03 1998-08-14 2021-03-24 1 49 313

This data set is long enough to approach the question of “Has there been a ‘measurable reduction in dissolved oxygen concentration’”?

First the assessment team divided the data into three depth bins:

  1. Upper (≤ 16.667 m)

  2. Mid (> 16.667 and ≤ 33.334 m)

  3. Bottom (> 33.334m)

Secondly, the data was divided into time periods, “IR period” representing the current integrated report window (≥ 2016-01-01) and historical data (< 2016-01-01).

newport_line_data <- newport_line_awqms %>%
  mutate(act_depth_height = as.numeric(act_depth_height)) %>%
  mutate(depth_position = case_when(act_depth_height <= 16.6667 ~ "Upper",
                                    act_depth_height <= 16.6667*2 ~ "Mid",
                                    act_depth_height <= 16.6667*3 ~ "Bottom")) %>%
  mutate(period = case_when(SampleStartDate < lubridate::ymd("2016-01-01") ~ "Historical",
                            SampleStartDate >= lubridate::ymd("2016-01-01") ~ "IR Period"))

The first analysis we did was to try to determine if the mean historical DO values are significantly different from the current IR data window means.

NH01

Site NH01 is the most inner site in the Newport line. This location has only 2 results in the “Bottom” depth bin, so we eliminated those from further analysis

newport_line_NH01 <- newport_line_data %>%
  filter(MLocID == "NH01",
         depth_position != "Bottom")

Next we want to get a sense of the normality of the data:

Due to the data’s bimodal distribution we chose to run a 2-sample Wilcoxon test to see if we can detect statistical significance between the two time periods. First we ran the test on the upper portion of the water column.

# Get 2 vectors of vector of results, depending on time period 
NH01_Upper_his <- newport_line_data %>%
  filter(MLocID == "NH01",
         depth_position == "Upper",
         period == 'Historical') %>%
  pull(Result_Numeric)


NH01_Upper_IR <- newport_line_data %>%
  filter(MLocID == "NH01",
         depth_position == "Upper",
         period == 'IR Period') %>%
  pull(Result_Numeric)



wilcox.test(NH01_Upper_his, NH01_Upper_IR,  alternative = "two.sided")
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  NH01_Upper_his and NH01_Upper_IR
## W = 2505092, p-value = 0.9485
## alternative hypothesis: true location shift is not equal to 0

This test, with a p-value of 0.948 suggests there is no significant difference between mean DO concentrations in the historical time period and the current IR window.

Next we looked at the Mid level data:

# Get 2 vectors of vector of results, depending on time period 
NH01_mid_his <- newport_line_data %>%
  filter(MLocID == "NH01",
         depth_position == "Mid",
         period == 'Historical') %>%
  pull(Result_Numeric)


NH01_mid_IR <- newport_line_data %>%
  filter(MLocID == "NH01",
         depth_position == "Mid",
         period == 'IR Period') %>%
  pull(Result_Numeric)



wilcox.test(NH01_mid_his, NH01_mid_IR,  alternative = "two.sided")
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  NH01_mid_his and NH01_mid_IR
## W = 511110, p-value = 0.7631
## alternative hypothesis: true location shift is not equal to 0

This test, with a p-value of 0.763 suggests there is no significant difference between mean DO concentrations in the historical time period and the current IR window.

A boxplot of this data confirms.

However, if we drill down into August specifically, we start to see some statistical differences in means.

NH01_Upper_his_aug <- newport_line_data %>%
  filter(MLocID == "NH01",
         depth_position == "Upper",
         period == 'Historical',
         month(SampleStartDate) == 8) %>%
  pull(Result_Numeric)


NH01_Upper_IR_aug <- newport_line_data %>%
  filter(MLocID == "NH01",
         depth_position == "Upper",
         period == 'IR Period',
         month(SampleStartDate) == 8) %>%
  pull(Result_Numeric)


wilcox.test(NH01_Upper_his_aug, NH01_Upper_IR_aug,  alternative = "two.sided")
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  NH01_Upper_his_aug and NH01_Upper_IR_aug
## W = 32307, p-value = 0.0000000000000008103
## alternative hypothesis: true location shift is not equal to 0

This test, with a p-value of 0.0000000000000008103 suggests there is is significant difference between mean DO concentrations in the historical time period and the current IR window in the month of august.

The integrated report window has lower mean do concentrations in both the Upper and Mid portions of the water column.

NH03

We repeated this same analysis for location NH03.

newport_line_NH03 <- newport_line_data %>%
  filter(MLocID == "NH03")

Again, we see a bi-modal distribution so we chose to run a 2-sample Wilcoxon test to see if we can detect statistical significance between the two time periods. First we ran the test on the upper portion of the water column.

# Get 2 vectors of vector of results, depending on time period 
NH03_Upper_his <- newport_line_data %>%
  filter(MLocID == "NH03",
         depth_position == "Upper",
         period == 'Historical') %>%
  pull(Result_Numeric)


NH03_Upper_IR <- newport_line_data %>%
  filter(MLocID == "NH03",
         depth_position == "Upper",
         period == 'IR Period') %>%
  pull(Result_Numeric)


wilcox.test(NH03_Upper_his, NH03_Upper_IR,  alternative = "two.sided")
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  NH03_Upper_his and NH03_Upper_IR
## W = 2481276, p-value = 0.8563
## alternative hypothesis: true location shift is not equal to 0

This test, with a p-value of 0.856 suggests there is no significant difference between mean DO concentrations in the historical time period and the current IR window.

Next we looked at the Mid level data:

# Get 2 vectors of vector of results, depending on time period 
NH03_mid_his <- newport_line_data %>%
  filter(MLocID == "NH03",
         depth_position == "Mid",
         period == 'Historical') %>%
  pull(Result_Numeric)


NH03_mid_IR <- newport_line_data %>%
  filter(MLocID == "NH03",
         depth_position == "Mid",
         period == 'IR Period') %>%
  pull(Result_Numeric)

wilcox.test(NH03_mid_his, NH03_mid_IR,  alternative = "two.sided")
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  NH03_mid_his and NH03_mid_IR
## W = 2864278, p-value = 0.03612
## alternative hypothesis: true location shift is not equal to 0

This test, with a p-value of 0.036 suggests there is a significant difference between mean DO concentrations in the historical time period and the current IR window.

And now for the bottom depth position:

# Get 2 vectors of vector of results, depending on time period 
NH03_bot_his <- newport_line_data %>%
  filter(MLocID == "NH03",
         depth_position == "Bottom",
         period == 'Historical') %>%
  pull(Result_Numeric)


NH03_bot_IR <- newport_line_data %>%
  filter(MLocID == "NH03",
         depth_position == "Bottom",
         period == 'IR Period') %>%
  pull(Result_Numeric)

wilcox.test(NH03_bot_his, NH03_bot_IR,  alternative = "two.sided")
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  NH03_bot_his and NH03_bot_IR
## W = 1354553, p-value = 0.03464
## alternative hypothesis: true location shift is not equal to 0

Again, this test, with a p-value of 0.035 suggests there is a significant difference between mean DO concentrations in the historical time period and the current IR window.

The boxplots show that while the means have a slight statistical significant difference, the relative distributions of the DO concentrations is roughly the same.

We look again at August specifically.

# Get 2 vectors of vector of results, depending on time period 
NH03_Upper_his_aug <- newport_line_data %>%
  filter(MLocID == "NH03",
         depth_position == "Upper",
         period == 'Historical', 
         month(SampleStartDate) == 8) %>%
  pull(Result_Numeric)


NH03_Upper_IR_aug <- newport_line_data %>%
  filter(MLocID == "NH03",
         depth_position == "Upper",
         period == 'IR Period', 
         month(SampleStartDate) == 8) %>%
  pull(Result_Numeric)


wilcox.test(NH03_Upper_his_aug, NH03_Upper_IR_aug,  alternative = "two.sided")
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  NH03_Upper_his_aug and NH03_Upper_IR_aug
## W = 32977, p-value = 0.0000001278
## alternative hypothesis: true location shift is not equal to 0

With a p-value of 0.0000001 we again see a statistical difference in the upper portion of the water column between historical data and the IR data window.

We see similar results in the mid and bottom level.

Mid Wilcoxon test results:

# Get 2 vectors of vector of results, depending on time period 
NH03_mid_his_aug <- newport_line_data %>%
  filter(MLocID == "NH03",
         depth_position == "Mid",
         period == 'Historical', 
         month(SampleStartDate) == 8) %>%
  pull(Result_Numeric)


NH03_mid_IR_aug <- newport_line_data %>%
  filter(MLocID == "NH03",
         depth_position == "Mid",
         period == 'IR Period', 
         month(SampleStartDate) == 8) %>%
  pull(Result_Numeric)


wilcox.test(NH03_mid_his_aug, NH03_mid_IR_aug,  alternative = "two.sided")
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  NH03_mid_his_aug and NH03_mid_IR_aug
## W = 41395, p-value = 0.00000000002786
## alternative hypothesis: true location shift is not equal to 0

An d the bottom Wilcoxon test results:

# Get 2 vectors of vector of results, depending on time period 
NH03_bottom_his_aug <- newport_line_data %>%
  filter(MLocID == "NH03",
         depth_position == "Bottom",
         period == 'Historical', 
         month(SampleStartDate) == 8) %>%
  pull(Result_Numeric)


NH03_bottom_IR_aug <- newport_line_data %>%
  filter(MLocID == "NH03",
         depth_position == "Bottom",
         period == 'IR Period', 
         month(SampleStartDate) == 8) %>%
  pull(Result_Numeric)


wilcox.test(NH03_bottom_his_aug, NH03_bottom_IR_aug,  alternative = "two.sided")
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  NH03_bottom_his_aug and NH03_bottom_IR_aug
## W = 19289, p-value = 0.0002316
## alternative hypothesis: true location shift is not equal to 0

The p-values for the mid and bottom are 0 and 0.0002316, respectively.

The box plot for NH03 in August is as follows.:

Next, we investigated all the months to see if there were differences in the months.

month_stats <- newport_line_data %>%
  filter(!(MLocID == "NH01" & depth_position == "Bottom")) %>%
  mutate(month =  month(SampleStartDate, label = TRUE)) %>%
  select(MLocID, depth_position, period, month, Result_Numeric) %>%
  group_by(MLocID, month, depth_position) %>%
  nest()%>%
  dplyr::summarise(wilcoxxon_pvalue = purrr::map(data, ~ .x %>%
                                                   summarise(wilcoxxon_pvalue = wilcox.test(Result_Numeric[period == "Historical"],
                                                                                            Result_Numeric[period == "IR Period"])[["p.value"]]))) %>%
  tidyr::unnest_wider(wilcoxxon_pvalue) %>%
  mutate(statistical_signifigance = case_when(wilcoxxon_pvalue < 0.05/12 ~ "Signifigant",
                                              wilcoxxon_pvalue >= 0.05/12 ~ "No statistical signifigance"))
p-value > 0.00416 is signifigant
MLocID month Mid Upper Bottom
NH01 Jan Signifigant Signifigant NA
NH01 Feb No statistical signifigance No statistical signifigance NA
NH01 Mar No statistical signifigance No statistical signifigance NA
NH01 Apr No statistical signifigance Signifigant NA
NH01 May Signifigant Signifigant NA
NH01 Jun No statistical signifigance No statistical signifigance NA
NH01 Jul No statistical signifigance No statistical signifigance NA
NH01 Aug Signifigant Signifigant NA
NH01 Sep Signifigant Signifigant NA
NH01 Oct No statistical signifigance No statistical signifigance NA
NH01 Nov No statistical signifigance No statistical signifigance NA
NH01 Dec No statistical signifigance No statistical signifigance NA
NH03 Jan No statistical signifigance Signifigant No statistical signifigance
NH03 Feb Signifigant Signifigant No statistical signifigance
NH03 Mar Signifigant Signifigant Signifigant
NH03 Apr No statistical signifigance Signifigant No statistical signifigance
NH03 May Signifigant No statistical signifigance No statistical signifigance
NH03 Jun No statistical signifigance No statistical signifigance Signifigant
NH03 Jul Signifigant No statistical signifigance Signifigant
NH03 Aug Signifigant Signifigant Signifigant
NH03 Sep Signifigant Signifigant Signifigant
NH03 Oct No statistical signifigance No statistical signifigance No statistical signifigance
NH03 Nov No statistical signifigance No statistical signifigance Signifigant
NH03 Dec No statistical signifigance Signifigant No statistical signifigance

To visualize all the months together:

and

Given the month to month differences in mean DO levels, the assessment teams concludes this assment unit should ramain categorized as category 3B- insufficient data; potential concern. Hopefully, our work group will be able to help DEQ understand the effects of marine DO, and acidification on fish and wildlife communities, as well as more fully understand long term trends in marine dissolved oxygen in territorial waters.